home *** CD-ROM | disk | FTP | other *** search
/ Cream of the Crop 1 / Cream of the Crop 1.iso / PROGRAM / DJLSR106.ARJ / INTEGER.CC < prev    next >
C/C++ Source or Header  |  1992-03-30  |  52KB  |  2,410 lines
  1. /* 
  2. Copyright (C) 1988 Free Software Foundation
  3.     written by Doug Lea (dl@rocky.oswego.edu)
  4.  
  5. This file is part of the GNU C++ Library.  This library is free
  6. software; you can redistribute it and/or modify it under the terms of
  7. the GNU Library General Public License as published by the Free
  8. Software Foundation; either version 2 of the License, or (at your
  9. option) any later version.  This library is distributed in the hope
  10. that it will be useful, but WITHOUT ANY WARRANTY; without even the
  11. implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR
  12. PURPOSE.  See the GNU Library General Public License for more details.
  13. You should have received a copy of the GNU Library General Public
  14. License along with this library; if not, write to the Free Software
  15. Foundation, 675 Mass Ave, Cambridge, MA 02139, USA.
  16. */
  17.  
  18. /*
  19.   Some of the following algorithms are very loosely based on those from 
  20.   MIT C-Scheme bignum.c, which is
  21.       Copyright (c) 1987 Massachusetts Institute of Technology
  22.  
  23.   with other guidance from Knuth, vol. 2
  24.  
  25.   Thanks to the creators of the algorithms.
  26. */
  27.  
  28. #ifdef __GNUG__
  29. #pragma implementation "Integer.h"
  30. #endif
  31. #include <Integer.h>
  32. #include <std.h>
  33. #include <ctype.h>
  34. #include <math.h>
  35. #include <values.h>
  36. #include <_Obstack.h>
  37. #include <AllocRing.h>
  38. #include <new.h>
  39. #include <builtin.h>
  40.  
  41. /*
  42.  Sizes of shifts for multiple-precision arithmetic.
  43.  These should not be changed unless Integer representation
  44.  as unsigned shorts is changed in the implementation files.
  45. */
  46.  
  47. #define I_SHIFT         SHORTBITS
  48. #define I_RADIX         ((unsigned long)(1L << I_SHIFT))
  49. #define I_MAXNUM        ((unsigned long)((I_RADIX - 1)))
  50. #define I_MINNUM        ((unsigned long)(I_RADIX >> 1))
  51. #define I_POSITIVE      1
  52. #define I_NEGATIVE      0
  53.  
  54. /* All routines assume SHORT_PER_LONG > 1 */
  55. #define SHORT_PER_LONG  ((unsigned)(((LONGBITS + SHORTBITS - 1) / SHORTBITS)))
  56. #define CHAR_PER_LONG   ((unsigned)(((LONGBITS + CHARBITS - 1) / CHARBITS)))
  57.  
  58. /*
  59.   minimum and maximum sizes for an IntRep
  60. */
  61.  
  62. #define MINIntRep_SIZE   16
  63. #define MAXIntRep_SIZE   I_MAXNUM
  64.  
  65. #ifndef MALLOC_MIN_OVERHEAD
  66. #define MALLOC_MIN_OVERHEAD 4
  67. #endif
  68.  
  69.  
  70. // utilities to extract and transfer bits 
  71.  
  72. // get low bits
  73.  
  74. inline static unsigned short extract(unsigned long x)
  75. {
  76.   return x & I_MAXNUM;
  77. }
  78.  
  79. // transfer high bits to low
  80.  
  81. inline static unsigned long down(unsigned long x)
  82. {
  83.   return (x >> I_SHIFT) & I_MAXNUM;
  84. }
  85.  
  86. // transfer low bits to high
  87.  
  88. inline static unsigned long up(unsigned long x)
  89. {
  90.   return x << I_SHIFT;
  91. }
  92.  
  93. // compare two equal-length reps
  94.  
  95. inline static int docmp(const unsigned short* x, const unsigned short* y, int l)
  96. {
  97.   int diff = 0;
  98.   const unsigned short* xs = &(x[l]);
  99.   const unsigned short* ys = &(y[l]);
  100.   while (l-- > 0 && (diff = (*--xs) - (*--ys)) == 0);
  101.   return diff;
  102. }
  103.  
  104. // figure out max length of result of +, -, etc.
  105.  
  106. inline static int calc_len(int len1, int len2, int pad)
  107. {
  108.   return (len1 >= len2)? len1 + pad : len2 + pad;
  109. }
  110.  
  111. // ensure len & sgn are correct
  112.  
  113. inline static void Icheck(IntRep* rep)
  114. {
  115.   int l = rep->len;
  116.   const unsigned short* p = &(rep->s[l]);
  117.   while (l > 0 && *--p == 0) --l;
  118.   if ((rep->len = l) == 0) rep->sgn = I_POSITIVE;
  119. }
  120.  
  121.  
  122. // zero out the end of a rep
  123.  
  124. inline static void Iclear_from(IntRep* rep, int p)
  125. {
  126.   unsigned short* cp = &(rep->s[p]);
  127.   const unsigned short* cf = &(rep->s[rep->len]);
  128.   while(cp < cf) *cp++ = 0;
  129. }
  130.  
  131. // copy parts of a rep
  132.  
  133. static inline void scpy(const unsigned short* src, unsigned short* dest,int nb)
  134. {
  135.   while (--nb >= 0) *dest++ = *src++;
  136. }
  137.  
  138. // make sure an argument is valid
  139.  
  140. static inline void nonnil(const IntRep* rep)
  141. {
  142.   if (rep == 0) 
  143.     (*lib_error_handler)("Integer", "operation on uninitialized Integer");
  144. }
  145.  
  146. // allocate a new Irep. Pad to something close to a power of two.
  147.  
  148. inline static IntRep* Inew(int newlen)
  149. {
  150.   unsigned int siz = sizeof(IntRep) + newlen * sizeof(short) + 
  151.     MALLOC_MIN_OVERHEAD;
  152.   unsigned int allocsiz = MINIntRep_SIZE;
  153.   while (allocsiz < siz) allocsiz <<= 1;  // find a power of 2
  154.   allocsiz -= MALLOC_MIN_OVERHEAD;
  155.   if (allocsiz >= MAXIntRep_SIZE * sizeof(short))
  156.     (*lib_error_handler)("Integer", "Requested length out of range");
  157.     
  158.   IntRep* rep = (IntRep *) new char[allocsiz];
  159.   rep->sz = (allocsiz - sizeof(IntRep) + sizeof(short)) / sizeof(short);
  160.   return rep;
  161. }
  162.  
  163. // allocate: use the bits in src if non-null, clear the rest
  164.  
  165. IntRep* Ialloc(IntRep* old, const unsigned short* src, int srclen, int newsgn,
  166.               int newlen)
  167. {
  168.   IntRep* rep;
  169.   if (old == 0 || newlen > old->sz)
  170.     rep = Inew(newlen);
  171.   else
  172.     rep = old;
  173.  
  174.   rep->len = newlen;
  175.   rep->sgn = newsgn;
  176.  
  177.   scpy(src, rep->s, srclen);
  178.   Iclear_from(rep, srclen);
  179.  
  180.   if (old != rep && old != 0) delete old;
  181.   return rep;
  182. }
  183.  
  184. // allocate and clear
  185.  
  186. IntRep* Icalloc(IntRep* old, int newlen)
  187. {
  188.   IntRep* rep;
  189.   if (old == 0 || newlen > old->sz)
  190.   {
  191.     if (old != 0) delete old;
  192.     rep = Inew(newlen);
  193.   }
  194.   else
  195.     rep = old;
  196.  
  197.   rep->len = newlen;
  198.   rep->sgn = I_POSITIVE;
  199.   Iclear_from(rep, 0);
  200.  
  201.   return rep;
  202. }
  203.  
  204. // reallocate
  205.  
  206. IntRep* Iresize(IntRep* old, int newlen)
  207. {
  208.   IntRep* rep;
  209.   unsigned short oldlen;
  210.   if (old == 0)
  211.   {
  212.     oldlen = 0;
  213.     rep = Inew(newlen);
  214.     rep->sgn = I_POSITIVE;
  215.   }
  216.   else 
  217.   {
  218.     oldlen = old->len;
  219.     if (newlen > old->sz)
  220.     {
  221.       rep = Inew(newlen);
  222.       scpy(old->s, rep->s, oldlen);
  223.       rep->sgn = old->sgn;
  224.       delete old;
  225.     }
  226.     else
  227.       rep = old;
  228.   }
  229.  
  230.   rep->len = newlen;
  231.   Iclear_from(rep, oldlen);
  232.  
  233.   return rep;
  234. }
  235.  
  236.  
  237. // same, for straight copy
  238.  
  239. IntRep* Icopy(IntRep* old, const IntRep* src)
  240. {
  241.   if (old == src) return old; 
  242.   IntRep* rep;
  243.   if (src == 0)
  244.   {
  245.     if (old == 0)
  246.       rep = Inew(0);
  247.     else
  248.     {
  249.       rep = old;
  250.       Iclear_from(rep, 0);
  251.     }
  252.     rep->len = 0;
  253.     rep->sgn = I_POSITIVE;
  254.   }
  255.   else 
  256.   {
  257.     int newlen = src->len;
  258.     if (old == 0 || newlen > old->sz)
  259.     {
  260.       if (old != 0) delete old;
  261.       rep = Inew(newlen);
  262.     }
  263.     else
  264.       rep = old;
  265.  
  266.     rep->len = newlen;
  267.     rep->sgn = src->sgn;
  268.  
  269.     scpy(src->s, rep->s, newlen);
  270.   }
  271.  
  272.   return rep;
  273. }
  274.  
  275. // allocate & copy space for a long
  276.  
  277. IntRep* Icopy_long(IntRep* old, long x)
  278. {
  279.   unsigned short src[SHORT_PER_LONG];
  280.   unsigned long u;
  281.   int newsgn;
  282.   if (newsgn = (x >= 0))
  283.     u = x;
  284.   else
  285.     u = -x;
  286.   
  287.   unsigned short srclen = 0;
  288.   while (u != 0)
  289.   {
  290.     src[srclen++] = extract(u);
  291.     u = down(u);
  292.   }
  293.  
  294.   IntRep* rep;
  295.   if (old == 0 || srclen > old->sz)
  296.   {
  297.     if (old != 0) delete old;
  298.     rep = Inew(srclen);
  299.   }
  300.   else
  301.     rep = old;
  302.  
  303.   rep->len = srclen;
  304.   rep->sgn = newsgn;
  305.  
  306.   scpy(src, rep->s, srclen);
  307.  
  308.   return rep;
  309.  
  310. }
  311.  
  312. // special case for zero -- it's worth it!
  313.  
  314. IntRep* Icopy_zero(IntRep* old)
  315. {
  316.   IntRep* rep;
  317.   if (old == 0)
  318.     rep = Inew(0);
  319.   else
  320.     rep = old;
  321.  
  322.   rep->len = 0;
  323.   rep->sgn = I_POSITIVE;
  324.  
  325.   return rep;
  326. }
  327.  
  328. // special case for 1 or -1
  329.  
  330. IntRep* Icopy_one(IntRep* old, int newsgn)
  331. {
  332.   IntRep* rep;
  333.   if (old == 0 || 1 > old->sz)
  334.   {
  335.     if (old != 0) delete old;
  336.     rep = Inew(1);
  337.   }
  338.   else
  339.     rep = old;
  340.  
  341.   rep->sgn = newsgn;
  342.   rep->len = 1;
  343.   rep->s[0] = 1;
  344.  
  345.   return rep;
  346. }
  347.  
  348. // convert to a legal two's complement long if possible
  349. // if too big, return most negative/positive value
  350.  
  351. long Itolong(const IntRep* rep)
  352. { 
  353.   if ((unsigned)(rep->len) > (unsigned)(SHORT_PER_LONG))
  354.     return (rep->sgn == I_POSITIVE) ? MAXLONG : MINLONG;
  355.   else if (rep->len == 0)
  356.     return 0;
  357.   else if ((unsigned)(rep->len) < (unsigned)(SHORT_PER_LONG))
  358.   {
  359.     unsigned long a = rep->s[rep->len-1];
  360.     if (SHORT_PER_LONG > 2) // normally optimized out
  361.     {
  362.       for (int i = rep->len - 2; i >= 0; --i)
  363.         a = up(a) | rep->s[i];
  364.     }
  365.     return (rep->sgn == I_POSITIVE)? a : -((long)a);
  366.   }
  367.   else 
  368.   {
  369.     unsigned long a = rep->s[SHORT_PER_LONG - 1];
  370.     if (a >= I_MINNUM)
  371.       return (rep->sgn == I_POSITIVE) ? MAXLONG : MINLONG;
  372.     else
  373.     {
  374.       a = up(a) | rep->s[SHORT_PER_LONG - 2];
  375.       if (SHORT_PER_LONG > 2)
  376.       {
  377.         for (int i = SHORT_PER_LONG - 3; i >= 0; --i)
  378.           a = up(a) | rep->s[i];
  379.       }
  380.       return (rep->sgn == I_POSITIVE)? a : -((long)a);
  381.     }
  382.   }
  383. }
  384.  
  385. // test whether op long() will work.
  386. // careful about asymmetry between MINLONG & MAXLONG
  387.  
  388. int Iislong(const IntRep* rep)
  389. {
  390.   unsigned int l = rep->len;
  391.   if (l < SHORT_PER_LONG)
  392.     return 1;
  393.   else if (l > SHORT_PER_LONG)
  394.     return 0;
  395.   else if ((unsigned)(rep->s[SHORT_PER_LONG - 1]) < (unsigned)(I_MINNUM))
  396.     return 1;
  397.   else if (rep->sgn == I_NEGATIVE && rep->s[SHORT_PER_LONG - 1] == I_MINNUM)
  398.   {
  399.     for (unsigned int i = 0; i < SHORT_PER_LONG - 1; ++i)
  400.       if (rep->s[i] != 0)
  401.         return 0;
  402.     return 1;
  403.   }
  404.   else
  405.     return 0;
  406. }
  407.  
  408. // convert to a double 
  409.  
  410. double Itodouble(const IntRep* rep)
  411. { 
  412.   double d = 0.0;
  413.   double bound = HUGE / 2.0;
  414.   for (int i = rep->len - 1; i >= 0; --i)
  415.   {
  416.     unsigned short a = I_RADIX >> 1;
  417.     while (a != 0)
  418.     {
  419.       if (d >= bound)
  420.         return (rep->sgn == I_NEGATIVE) ? -HUGE : HUGE;
  421.       d *= 2.0;
  422.       if (rep->s[i] & a)
  423.         d += 1.0;
  424.       a >>= 1;
  425.     }
  426.   }
  427.   if (rep->sgn == I_NEGATIVE)
  428.     return -d;
  429.   else
  430.     return d;
  431. }
  432.  
  433. // see whether op double() will work-
  434. // have to actually try it in order to find out
  435. // since otherwise might trigger fp exception
  436.  
  437. int Iisdouble(const IntRep* rep)
  438. {
  439.   double d = 0.0;
  440.   double bound = HUGE / 2.0;
  441.   for (int i = rep->len - 1; i >= 0; --i)
  442.   {
  443.     unsigned short a = I_RADIX >> 1;
  444.     while (a != 0)
  445.     {
  446.       if (d > bound || (d == bound && (i > 0 || (rep->s[i] & a))))
  447.         return 0;
  448.       d *= 2.0;
  449.       if (rep->s[i] & a)
  450.         d += 1.0;
  451.       a >>= 1;
  452.     }
  453.   }
  454.   return 1;
  455. }
  456.  
  457. // real division of num / den
  458.  
  459. double ratio(const Integer& num, const Integer& den)
  460. {
  461.   Integer q, r;
  462.   divide(num, den, q, r);
  463.   double d1 = double(q);
  464.  
  465.   if (d1 == HUGE || d1 == -HUGE || sign(r) == 0)
  466.     return d1;
  467.   else      // use as much precision as available for fractional part
  468.   {
  469.     double  d2 = 0.0;
  470.     double  d3 = 0.0; 
  471.     int cont = 1;
  472.     for (int i = den.rep->len - 1; i >= 0 && cont; --i)
  473.     {
  474.       unsigned short a = I_RADIX >> 1;
  475.       while (a != 0)
  476.       {
  477.         if (d2 + 1.0 == d2) // out of precision when we get here
  478.         {
  479.           cont = 0;
  480.           break;
  481.         }
  482.  
  483.         d2 *= 2.0;
  484.         if (den.rep->s[i] & a)
  485.           d2 += 1.0;
  486.  
  487.         if (i < r.rep->len)
  488.         {
  489.           d3 *= 2.0;
  490.           if (r.rep->s[i] & a)
  491.             d3 += 1.0;
  492.         }
  493.  
  494.         a >>= 1;
  495.       }
  496.     }
  497.  
  498.     if (sign(r) < 0)
  499.       d3 = -d3;
  500.     return d1 + d3 / d2;
  501.   }
  502. }
  503.  
  504. // comparison functions
  505.   
  506. int compare(const IntRep* x, const IntRep* y)
  507. {
  508.   int diff  = x->sgn - y->sgn;
  509.   if (diff == 0)
  510.   {
  511.     diff = x->len - y->len;
  512.     if (diff == 0)
  513.       diff = docmp(x->s, y->s, x->len);
  514.     if (x->sgn == I_NEGATIVE)
  515.       diff = -diff;
  516.   }
  517.   return diff;
  518. }
  519.  
  520. int ucompare(const IntRep* x, const IntRep* y)
  521. {
  522.   int diff = x->len - y->len;
  523.   if (diff == 0)
  524.   {
  525.     int l = x->len;
  526.     const unsigned short* xs = &(x->s[l]);
  527.     const unsigned short* ys = &(y->s[l]);
  528.     while (l-- > 0 && (diff = (*--xs) - (*--ys)) == 0);
  529.   }
  530.   return diff;
  531. }
  532.  
  533. int compare(const IntRep* x, long  y)
  534. {
  535.   int xl = x->len;
  536.   int xsgn = x->sgn;
  537.   if (y == 0)
  538.   {
  539.     if (xl == 0)
  540.       return 0;
  541.     else if (xsgn == I_NEGATIVE)
  542.       return -1;
  543.     else
  544.       return 1;
  545.   }
  546.   else
  547.   {
  548.     int ysgn = y >= 0;
  549.     unsigned long uy = (ysgn)? y : -y;
  550.     int diff = xsgn - ysgn;
  551.     if (diff == 0)
  552.     {
  553.       diff = xl - SHORT_PER_LONG;
  554.       if (diff <= 0)
  555.       {
  556.         unsigned short tmp[SHORT_PER_LONG];
  557.         int yl = 0;
  558.         while (uy != 0)
  559.         {
  560.           tmp[yl++] = extract(uy);
  561.           uy = down(uy);
  562.         }
  563.         diff = xl - yl;
  564.         if (diff == 0)
  565.           diff = docmp(x->s, tmp, xl);
  566.       }
  567.     }
  568.     if (xsgn == I_NEGATIVE)
  569.       diff = -diff;
  570.     return diff;
  571.   }
  572. }
  573.  
  574. int ucompare(const IntRep* x, long  y)
  575. {
  576.   int xl = x->len;
  577.   if (y == 0)
  578.     return xl;
  579.   else
  580.   {
  581.     unsigned long uy = (y >= 0)? y : -y;
  582.     int diff = xl - SHORT_PER_LONG;
  583.     if (diff <= 0)
  584.     {
  585.       unsigned short tmp[SHORT_PER_LONG];
  586.       int yl = 0;
  587.       while (uy != 0)
  588.       {
  589.         tmp[yl++] = extract(uy);
  590.         uy = down(uy);
  591.       }
  592.       diff = xl - yl;
  593.       if (diff == 0)
  594.         diff = docmp(x->s, tmp, xl);
  595.     }
  596.     return diff;
  597.   }
  598. }
  599.  
  600.  
  601.  
  602. // arithmetic functions
  603.  
  604. IntRep* add(const IntRep* x, int negatex, 
  605.             const IntRep* y, int negatey, IntRep* r)
  606. {
  607.   nonnil(x);
  608.   nonnil(y);
  609.  
  610.   int xl = x->len;
  611.   int yl = y->len;
  612.  
  613.   int xsgn = (negatex && xl != 0) ? !x->sgn : x->sgn;
  614.   int ysgn = (negatey && yl != 0) ? !y->sgn : y->sgn;
  615.  
  616.   int xrsame = x == r;
  617.   int yrsame = y == r;
  618.  
  619.   if (yl == 0)
  620.     r = Ialloc(r, x->s, xl, xsgn, xl);
  621.   else if (xl == 0)
  622.     r = Ialloc(r, y->s, yl, ysgn, yl);
  623.   else if (xsgn == ysgn)
  624.   {
  625.     if (xrsame || yrsame)
  626.       r = Iresize(r, calc_len(xl, yl, 1));
  627.     else
  628.       r = Icalloc(r, calc_len(xl, yl, 1));
  629.     r->sgn = xsgn;
  630.     unsigned short* rs = r->s;
  631.     const unsigned short* as;
  632.     const unsigned short* bs;
  633.     const unsigned short* topa;
  634.     const unsigned short* topb;
  635.     if (xl >= yl)
  636.     {
  637.       as =  (xrsame)? r->s : x->s;
  638.       topa = &(as[xl]);
  639.       bs =  (yrsame)? r->s : y->s;
  640.       topb = &(bs[yl]);
  641.     }
  642.     else
  643.     {
  644.       bs =  (xrsame)? r->s : x->s;
  645.       topb = &(bs[xl]);
  646.       as =  (yrsame)? r->s : y->s;
  647.       topa = &(as[yl]);
  648.     }
  649.     unsigned long sum = 0;
  650.     while (bs < topb)
  651.     {
  652.       sum += (unsigned long)(*as++) + (unsigned long)(*bs++);
  653.       *rs++ = extract(sum);
  654.       sum = down(sum);
  655.     }
  656.     while (sum != 0 && as < topa)
  657.     {
  658.       sum += (unsigned long)(*as++);
  659.       *rs++ = extract(sum);
  660.       sum = down(sum);
  661.     }
  662.     if (sum != 0)
  663.       *rs = extract(sum);
  664.     else if (rs != as)
  665.       while (as < topa)
  666.         *rs++ = *as++;
  667.   }
  668.   else
  669.   {
  670.     int comp = ucompare(x, y);
  671.     if (comp == 0)
  672.       r = Icopy_zero(r);
  673.     else
  674.     {
  675.       if (xrsame || yrsame)
  676.         r = Iresize(r, calc_len(xl, yl, 0));
  677.       else
  678.         r = Icalloc(r, calc_len(xl, yl, 0));
  679.       unsigned short* rs = r->s;
  680.       const unsigned short* as;
  681.       const unsigned short* bs;
  682.       const unsigned short* topa;
  683.       const unsigned short* topb;
  684.       if (comp > 0)
  685.       {
  686.         as =  (xrsame)? r->s : x->s;
  687.         topa = &(as[xl]);
  688.         bs =  (yrsame)? r->s : y->s;
  689.         topb = &(bs[yl]);
  690.         r->sgn = xsgn;
  691.       }
  692.       else
  693.       {
  694.         bs =  (xrsame)? r->s : x->s;
  695.         topb = &(bs[xl]);
  696.         as =  (yrsame)? r->s : y->s;
  697.         topa = &(as[yl]);
  698.         r->sgn = ysgn;
  699.       }
  700.       unsigned long hi = 1;
  701.       while (bs < topb)
  702.       {
  703.         hi += (unsigned long)(*as++) + I_MAXNUM - (unsigned long)(*bs++);
  704.         *rs++ = extract(hi);
  705.         hi = down(hi);
  706.       }
  707.       while (hi == 0 && as < topa)
  708.       {
  709.         hi = (unsigned long)(*as++) + I_MAXNUM;
  710.         *rs++ = extract(hi);
  711.         hi = down(hi);
  712.       }
  713.       if (rs != as)
  714.         while (as < topa)
  715.           *rs++ = *as++;
  716.     }
  717.   }
  718.   Icheck(r);
  719.   return r;
  720. }
  721.  
  722.  
  723. IntRep* add(const IntRep* x, int negatex, long y, IntRep* r)
  724. {
  725.   nonnil(x);
  726.   int xl = x->len;
  727.   int xsgn = (negatex && xl != 0) ? !x->sgn : x->sgn;
  728.   int xrsame = x == r;
  729.  
  730.   int ysgn = (y >= 0);
  731.   unsigned long uy = (ysgn)? y : -y;
  732.  
  733.   if (y == 0)
  734.     r = Ialloc(r, x->s, xl, xsgn, xl);
  735.   else if (xl == 0)
  736.     r = Icopy_long(r, y);
  737.   else if (xsgn == ysgn)
  738.   {
  739.     if (xrsame)
  740.       r = Iresize(r, calc_len(xl, SHORT_PER_LONG, 1));
  741.     else
  742.       r = Icalloc(r, calc_len(xl, SHORT_PER_LONG, 1));
  743.     r->sgn = xsgn;
  744.     unsigned short* rs = r->s;
  745.     const unsigned short* as =  (xrsame)? r->s : x->s;
  746.     const unsigned short* topa = &(as[xl]);
  747.     unsigned long sum = 0;
  748.     while (as < topa && uy != 0)
  749.     {
  750.       unsigned long u = extract(uy);
  751.       uy = down(uy);
  752.       sum += (unsigned long)(*as++) + u;
  753.       *rs++ = extract(sum);
  754.       sum = down(sum);
  755.     }
  756.     while (sum != 0 && as < topa)
  757.     {
  758.       sum += (unsigned long)(*as++);
  759.       *rs++ = extract(sum);
  760.       sum = down(sum);
  761.     }
  762.     if (sum != 0)
  763.       *rs = extract(sum);
  764.     else if (rs != as)
  765.       while (as < topa)
  766.         *rs++ = *as++;
  767.   }
  768.   else
  769.   {
  770.     unsigned short tmp[SHORT_PER_LONG];
  771.     int yl = 0;
  772.     while (uy != 0)
  773.     {
  774.       tmp[yl++] = extract(uy);
  775.       uy = down(uy);
  776.     }
  777.     int comp = xl - yl;
  778.     if (comp == 0)
  779.       comp = docmp(x->s, tmp, yl);
  780.     if (comp == 0)
  781.       r = Icopy_zero(r);
  782.     else
  783.     {
  784.       if (xrsame)
  785.         r = Iresize(r, calc_len(xl, yl, 0));
  786.       else
  787.         r = Icalloc(r, calc_len(xl, yl, 0));
  788.       unsigned short* rs = r->s;
  789.       const unsigned short* as;
  790.       const unsigned short* bs;
  791.       const unsigned short* topa;
  792.       const unsigned short* topb;
  793.       if (comp > 0)
  794.       {
  795.         as =  (xrsame)? r->s : x->s;
  796.         topa = &(as[xl]);
  797.         bs =  tmp;
  798.         topb = &(bs[yl]);
  799.         r->sgn = xsgn;
  800.       }
  801.       else
  802.       {
  803.         bs =  (xrsame)? r->s : x->s;
  804.         topb = &(bs[xl]);
  805.         as =  tmp;
  806.         topa = &(as[yl]);
  807.         r->sgn = ysgn;
  808.       }
  809.       unsigned long hi = 1;
  810.       while (bs < topb)
  811.       {
  812.         hi += (unsigned long)(*as++) + I_MAXNUM - (unsigned long)(*bs++);
  813.         *rs++ = extract(hi);
  814.         hi = down(hi);
  815.       }
  816.       while (hi == 0 && as < topa)
  817.       {
  818.         hi = (unsigned long)(*as++) + I_MAXNUM;
  819.         *rs++ = extract(hi);
  820.         hi = down(hi);
  821.       }
  822.       if (rs != as)
  823.         while (as < topa)
  824.           *rs++ = *as++;
  825.     }
  826.   }
  827.   Icheck(r);
  828.   return r;
  829. }
  830.  
  831.  
  832. IntRep* multiply(const IntRep* x, const IntRep* y, IntRep* r)
  833. {
  834.   nonnil(x);
  835.   nonnil(y);
  836.   int xl = x->len;
  837.   int yl = y->len;
  838.   int rl = xl + yl;
  839.   int rsgn = x->sgn == y->sgn;
  840.   int xrsame = x == r;
  841.   int yrsame = y == r;
  842.   int xysame = x == y;
  843.   
  844.   if (xl == 0 || yl == 0)
  845.     r = Icopy_zero(r);
  846.   else if (xl == 1 && x->s[0] == 1)
  847.     r = Icopy(r, y);
  848.   else if (yl == 1 && y->s[0] == 1)
  849.     r = Icopy(r, x);
  850.   else if (!(xysame && xrsame))
  851.   {
  852.     if (xrsame || yrsame)
  853.       r = Iresize(r, rl);
  854.     else
  855.       r = Icalloc(r, rl);
  856.     unsigned short* rs = r->s;
  857.     unsigned short* topr = &(rs[rl]);
  858.  
  859.     // use best inner/outer loop params given constraints
  860.     unsigned short* currentr;
  861.     const unsigned short* bota;
  862.     const unsigned short* as;
  863.     const unsigned short* botb;
  864.     const unsigned short* topb;
  865.     if (xrsame)                 
  866.     { 
  867.       currentr = &(rs[xl-1]);
  868.       bota = rs;
  869.       as = currentr;
  870.       botb = y->s;
  871.       topb = &(botb[yl]);
  872.     }
  873.     else if (yrsame)
  874.     {
  875.       currentr = &(rs[yl-1]);
  876.       bota = rs;
  877.       as = currentr;
  878.       botb = x->s;
  879.       topb = &(botb[xl]);
  880.     }
  881.     else if (xl <= yl)
  882.     {
  883.       currentr = &(rs[xl-1]);
  884.       bota = x->s;
  885.       as = &(bota[xl-1]);
  886.       botb = y->s;
  887.       topb = &(botb[yl]);
  888.     }
  889.     else
  890.     {
  891.       currentr = &(rs[yl-1]);
  892.       bota = y->s;
  893.       as = &(bota[yl-1]);
  894.       botb = x->s;
  895.       topb = &(botb[xl]);
  896.     }
  897.  
  898.     while (as >= bota)
  899.     {
  900.       unsigned long ai = (unsigned long)(*as--);
  901.       unsigned short* rs = currentr--;
  902.       *rs = 0;
  903.       if (ai != 0)
  904.       {
  905.         unsigned long sum = 0;
  906.         const unsigned short* bs = botb;
  907.         while (bs < topb)
  908.         {
  909.           sum += ai * (unsigned long)(*bs++) + (unsigned long)(*rs);
  910.           *rs++ = extract(sum);
  911.           sum = down(sum);
  912.         }
  913.         while (sum != 0 && rs < topr)
  914.         {
  915.           sum += (unsigned long)(*rs);
  916.           *rs++ = extract(sum);
  917.           sum = down(sum);
  918.         }
  919.       }
  920.     }
  921.   }
  922.   else                          // x, y, and r same; compute over diagonals
  923.   {
  924.     r = Iresize(r, rl);
  925.     unsigned short* botr = r->s;
  926.     unsigned short* topr = &(botr[rl]);
  927.     unsigned short* rs =   &(botr[rl - 2]);
  928.  
  929.     const unsigned short* bota = (xrsame)? botr : x->s;
  930.     const unsigned short* loa =  &(bota[xl - 1]);
  931.     const unsigned short* hia =  loa;
  932.  
  933.     for (; rs >= botr; --rs)
  934.     {
  935.       const unsigned short* h = hia;
  936.       const unsigned short* l = loa;
  937.       unsigned long prod = (unsigned long)(*h) * (unsigned long)(*l);
  938.       *rs = 0;
  939.  
  940.       for(;;)
  941.       {
  942.         unsigned short* rt = rs;
  943.         unsigned long sum = prod + (unsigned long)(*rt);
  944.         *rt++ = extract(sum);
  945.         sum = down(sum);
  946.         while (sum != 0 && rt < topr)
  947.         {
  948.           sum += (unsigned long)(*rt);
  949.           *rt++ = extract(sum);
  950.           sum = down(sum);
  951.         }
  952.         if (h > l)
  953.         {
  954.           rt = rs;
  955.           sum = prod + (unsigned long)(*rt);
  956.           *rt++ = extract(sum);
  957.           sum = down(sum);
  958.           while (sum != 0 && rt < topr)
  959.           {
  960.             sum += (unsigned long)(*rt);
  961.             *rt++ = extract(sum);
  962.             sum = down(sum);
  963.           }
  964.           if (--h >= ++l)
  965.             prod = (unsigned long)(*h) * (unsigned long)(*l);
  966.           else
  967.             break;
  968.         }
  969.         else
  970.           break;
  971.       }
  972.       if (loa > bota)
  973.         --loa;
  974.       else
  975.         --hia;
  976.     }
  977.   }
  978.   r->sgn = rsgn;
  979.   Icheck(r);
  980.   return r;
  981. }
  982.  
  983.  
  984. IntRep* multiply(const IntRep* x, long y, IntRep* r)
  985. {
  986.   nonnil(x);
  987.   int xl = x->len;
  988.     
  989.   if (xl == 0 || y == 0)
  990.     r = Icopy_zero(r);
  991.   else if (y == 1)
  992.     r = Icopy(r, x);
  993.   else
  994.   {
  995.     int ysgn = y >= 0;
  996.     int rsgn = x->sgn == ysgn;
  997.     unsigned long uy = (ysgn)? y : -y;
  998.     unsigned short tmp[SHORT_PER_LONG];
  999.     int yl = 0;
  1000.     while (uy != 0)
  1001.     {
  1002.       tmp[yl++] = extract(uy);
  1003.       uy = down(uy);
  1004.     }
  1005.  
  1006.     int rl = xl + yl;
  1007.     int xrsame = x == r;
  1008.     if (xrsame)
  1009.       r = Iresize(r, rl);
  1010.     else
  1011.       r = Icalloc(r, rl);
  1012.  
  1013.     unsigned short* rs = r->s;
  1014.     unsigned short* topr = &(rs[rl]);
  1015.     unsigned short* currentr;
  1016.     const unsigned short* bota;
  1017.     const unsigned short* as;
  1018.     const unsigned short* botb;
  1019.     const unsigned short* topb;
  1020.  
  1021.     if (xrsame)
  1022.     { 
  1023.       currentr = &(rs[xl-1]);
  1024.       bota = rs;
  1025.       as = currentr;
  1026.       botb = tmp;
  1027.       topb = &(botb[yl]);
  1028.     }
  1029.     else if (xl <= yl)
  1030.     {
  1031.       currentr = &(rs[xl-1]);
  1032.       bota = x->s;
  1033.       as = &(bota[xl-1]);
  1034.       botb = tmp;
  1035.       topb = &(botb[yl]);
  1036.     }
  1037.     else
  1038.     {
  1039.       currentr = &(rs[yl-1]);
  1040.       bota = tmp;
  1041.       as = &(bota[yl-1]);
  1042.       botb = x->s;
  1043.       topb = &(botb[xl]);
  1044.     }
  1045.  
  1046.     while (as >= bota)
  1047.     {
  1048.       unsigned long ai = (unsigned long)(*as--);
  1049.       unsigned short* rs = currentr--;
  1050.       *rs = 0;
  1051.       if (ai != 0)
  1052.       {
  1053.         unsigned long sum = 0;
  1054.         const unsigned short* bs = botb;
  1055.         while (bs < topb)
  1056.         {
  1057.           sum += ai * (unsigned long)(*bs++) + (unsigned long)(*rs);
  1058.           *rs++ = extract(sum);
  1059.           sum = down(sum);
  1060.         }
  1061.         while (sum != 0 && rs < topr)
  1062.         {
  1063.           sum += (unsigned long)(*rs);
  1064.           *rs++ = extract(sum);
  1065.           sum = down(sum);
  1066.         }
  1067.       }
  1068.     }
  1069.     r->sgn = rsgn;
  1070.   }
  1071.   Icheck(r);
  1072.   return r;
  1073. }
  1074.  
  1075.  
  1076. // main division routine
  1077.  
  1078. static void do_divide(unsigned short* rs,
  1079.                       const unsigned short* ys, int yl,
  1080.                       unsigned short* qs, int ql)
  1081. {
  1082.   const unsigned short* topy = &(ys[yl]);
  1083.   unsigned short d1 = ys[yl - 1];
  1084.   unsigned short d2 = ys[yl - 2];
  1085.  
  1086.   int l = ql - 1;
  1087.   int i = l + yl;
  1088.   
  1089.   for (; l >= 0; --l, --i)
  1090.   {
  1091.     unsigned short qhat;       // guess q
  1092.     if (d1 == rs[i])
  1093.       qhat = I_MAXNUM;
  1094.     else
  1095.     {
  1096.       unsigned long lr = up((unsigned long)rs[i]) | rs[i-1];
  1097.       qhat = lr / d1;
  1098.     }
  1099.  
  1100.     for(;;)     // adjust q, use docmp to avoid overflow problems
  1101.     {
  1102.       unsigned short ts[3];
  1103.       unsigned long prod = (unsigned long)d2 * (unsigned long)qhat;
  1104.       ts[0] = extract(prod);
  1105.       prod = down(prod) + (unsigned long)d1 * (unsigned long)qhat;
  1106.       ts[1] = extract(prod);
  1107.       ts[2] = extract(down(prod));
  1108.       if (docmp(ts, &(rs[i-2]), 3) > 0)
  1109.         --qhat;
  1110.       else
  1111.         break;
  1112.     };
  1113.     
  1114.     // multiply & subtract
  1115.     
  1116.     const unsigned short* yt = ys;
  1117.     unsigned short* rt = &(rs[l]);
  1118.     unsigned long prod = 0;
  1119.     unsigned long hi = 1;
  1120.     while (yt < topy)
  1121.     {
  1122.       prod = (unsigned long)qhat * (unsigned long)(*yt++) + down(prod);
  1123.       hi += (unsigned long)(*rt) + I_MAXNUM - (unsigned long)(extract(prod));
  1124.       *rt++ = extract(hi);
  1125.       hi = down(hi);
  1126.     }
  1127.     hi += (unsigned long)(*rt) + I_MAXNUM - (unsigned long)(down(prod));
  1128.     *rt = extract(hi);
  1129.     hi = down(hi);
  1130.     
  1131.     // off-by-one, add back
  1132.     
  1133.     if (hi == 0)
  1134.     {
  1135.       --qhat;
  1136.       yt = ys;
  1137.       rt = &(rs[l]);
  1138.       hi = 0;
  1139.       while (yt < topy)
  1140.       {
  1141.         hi = (unsigned long)(*rt) + (unsigned long)(*yt++) + down(hi);
  1142.         *rt++ = extract(hi);
  1143.       }
  1144.       *rt = 0;
  1145.     }
  1146.     if (qs != 0)
  1147.       qs[l] = qhat;
  1148.   }
  1149. }
  1150.  
  1151. // divide by single digit, return remainder
  1152. // if q != 0, then keep the result in q, else just compute rem
  1153.  
  1154. static int unscale(const unsigned short* x, int xl, unsigned short y,
  1155.                    unsigned short* q)
  1156. {
  1157.   if (xl == 0 || y == 1)
  1158.     return 0;
  1159.   else if (q != 0)
  1160.   {
  1161.     unsigned short* botq = q;
  1162.     unsigned short* qs = &(botq[xl - 1]);
  1163.     const unsigned short* xs = &(x[xl - 1]);
  1164.     unsigned long rem = 0;
  1165.     while (qs >= botq)
  1166.     {
  1167.       rem = up(rem) | *xs--;
  1168.       unsigned long u = rem / y;
  1169.       *qs-- = extract(u);
  1170.       rem -= u * y;
  1171.     }
  1172.     int r = extract(rem);
  1173.     return r;
  1174.   }
  1175.   else                          // same loop, a bit faster if just need rem
  1176.   {
  1177.     const unsigned short* botx = x;
  1178.     const unsigned short* xs = &(botx[xl - 1]);
  1179.     unsigned long rem = 0;
  1180.     while (xs >= botx)
  1181.     {
  1182.       rem = up(rem) | *xs--;
  1183.       unsigned long u = rem / y;
  1184.       rem -= u * y;
  1185.     }
  1186.     int r = extract(rem);
  1187.     return r;
  1188.   }
  1189. }
  1190.  
  1191.  
  1192. IntRep* div(const IntRep* x, const IntRep* y, IntRep* q)
  1193. {
  1194.   nonnil(x);
  1195.   nonnil(y);
  1196.   int xl = x->len;
  1197.   int yl = y->len;
  1198.   if (yl == 0) (*lib_error_handler)("Integer", "attempted division by zero");
  1199.  
  1200.   int comp = ucompare(x, y);
  1201.   int xsgn = x->sgn;
  1202.   int ysgn = y->sgn;
  1203.  
  1204.   int samesign = xsgn == ysgn;
  1205.  
  1206.   if (comp < 0)
  1207.     q = Icopy_zero(q);
  1208.   else if (comp == 0)
  1209.     q = Icopy_one(q, samesign);
  1210.   else if (yl == 1)
  1211.   {
  1212.     q = Icopy(q, x);
  1213.     unscale(q->s, q->len, y->s[0], q->s);
  1214.   }
  1215.   else
  1216.   {
  1217.     IntRep* yy = 0;
  1218.     IntRep* r  = 0;
  1219.     unsigned short prescale = (I_RADIX / (1 + y->s[yl - 1]));
  1220.     if (prescale != 1 || y == q)
  1221.     {
  1222.       yy = multiply(y, ((long)prescale & I_MAXNUM), yy);
  1223.       r = multiply(x, ((long)prescale & I_MAXNUM), r);
  1224.     }
  1225.     else
  1226.     {
  1227.       yy = (IntRep*)y;
  1228.       r = Icalloc(r, xl + 1);
  1229.       scpy(x->s, r->s, xl);
  1230.     }
  1231.  
  1232.     int ql = xl - yl + 1;
  1233.       
  1234.     q = Icalloc(q, ql);
  1235.     do_divide(r->s, yy->s, yl, q->s, ql);
  1236.  
  1237.     if (yy != y) delete yy;
  1238.     delete r;
  1239.   }
  1240.   q->sgn = samesign;
  1241.   Icheck(q);
  1242.   return q;
  1243. }
  1244.  
  1245. IntRep* div(const IntRep* x, long y, IntRep* q)
  1246. {
  1247.   nonnil(x);
  1248.   int xl = x->len;
  1249.   if (y == 0) (*lib_error_handler)("Integer", "attempted division by zero");
  1250.  
  1251.   unsigned short ys[SHORT_PER_LONG];
  1252.   unsigned long u;
  1253.   int ysgn = y >= 0;
  1254.   if (ysgn)
  1255.     u = y;
  1256.   else
  1257.     u = -y;
  1258.   int yl = 0;
  1259.   while (u != 0)
  1260.   {
  1261.     ys[yl++] = extract(u);
  1262.     u = down(u);
  1263.   }
  1264.  
  1265.   int comp = xl - yl;
  1266.   if (comp == 0) comp = docmp(x->s, ys, xl);
  1267.  
  1268.   int xsgn = x->sgn;
  1269.   int samesign = xsgn == ysgn;
  1270.  
  1271.   if (comp < 0)
  1272.     q = Icopy_zero(q);
  1273.   else if (comp == 0)
  1274.   {
  1275.     q = Icopy_one(q, samesign);
  1276.   }
  1277.   else if (yl == 1)
  1278.   {
  1279.     q = Icopy(q, x);
  1280.     unscale(q->s, q->len, ys[0], q->s);
  1281.   }
  1282.   else
  1283.   {
  1284.     IntRep* r  = 0;
  1285.     unsigned short prescale = (I_RADIX / (1 + ys[yl - 1]));
  1286.     if (prescale != 1)
  1287.     {
  1288.       unsigned long prod = (unsigned long)prescale * (unsigned long)ys[0];
  1289.       ys[0] = extract(prod);
  1290.       prod = down(prod) + (unsigned long)prescale * (unsigned long)ys[1];
  1291.       ys[1] = extract(prod);
  1292.       r = multiply(x, ((long)prescale & I_MAXNUM), r);
  1293.     }
  1294.     else
  1295.     {
  1296.       r = Icalloc(r, xl + 1);
  1297.       scpy(x->s, r->s, xl);
  1298.     }
  1299.  
  1300.     int ql = xl - yl + 1;
  1301.       
  1302.     q = Icalloc(q, ql);
  1303.     do_divide(r->s, ys, yl, q->s, ql);
  1304.  
  1305.     delete r;
  1306.   }
  1307.   q->sgn = samesign;
  1308.   Icheck(q);
  1309.   return q;
  1310. }
  1311.  
  1312.  
  1313. void divide(const Integer& Ix, long y, Integer& Iq, long& rem)
  1314. {
  1315.   const IntRep* x = Ix.rep;
  1316.   nonnil(x);
  1317.   IntRep* q = Iq.rep;
  1318.   int xl = x->len;
  1319.   if (y == 0) (*lib_error_handler)("Integer", "attempted division by zero");
  1320.  
  1321.   unsigned short ys[SHORT_PER_LONG];
  1322.   unsigned long u;
  1323.   int ysgn = y >= 0;
  1324.   if (ysgn)
  1325.     u = y;
  1326.   else
  1327.     u = -y;
  1328.   int yl = 0;
  1329.   while (u != 0)
  1330.   {
  1331.     ys[yl++] = extract(u);
  1332.     u = down(u);
  1333.   }
  1334.  
  1335.   int comp = xl - yl;
  1336.   if (comp == 0) comp = docmp(x->s, ys, xl);
  1337.  
  1338.   int xsgn = x->sgn;
  1339.   int samesign = xsgn == ysgn;
  1340.  
  1341.   if (comp < 0)
  1342.   {
  1343.     rem = Itolong(x);
  1344.     q = Icopy_zero(q);
  1345.   }
  1346.   else if (comp == 0)
  1347.   {
  1348.     q = Icopy_one(q, samesign);
  1349.     rem = 0;
  1350.   }
  1351.   else if (yl == 1)
  1352.   {
  1353.     q = Icopy(q, x);
  1354.     rem = unscale(q->s, q->len, ys[0], q->s);
  1355.   }
  1356.   else
  1357.   {
  1358.     IntRep* r  = 0;
  1359.     unsigned short prescale = (I_RADIX / (1 + ys[yl - 1]));
  1360.     if (prescale != 1)
  1361.     {
  1362.       unsigned long prod = (unsigned long)prescale * (unsigned long)ys[0];
  1363.       ys[0] = extract(prod);
  1364.       prod = down(prod) + (unsigned long)prescale * (unsigned long)ys[1];
  1365.       ys[1] = extract(prod);
  1366.       r = multiply(x, ((long)prescale & I_MAXNUM), r);
  1367.     }
  1368.     else
  1369.     {
  1370.       r = Icalloc(r, xl + 1);
  1371.       scpy(x->s, r->s, xl);
  1372.     }
  1373.  
  1374.     int ql = xl - yl + 1;
  1375.       
  1376.     q = Icalloc(q, ql);
  1377.     
  1378.     do_divide(r->s, ys, yl, q->s, ql);
  1379.  
  1380.     if (prescale != 1)
  1381.     {
  1382.       Icheck(r);
  1383.       unscale(r->s, r->len, prescale, r->s);
  1384.     }
  1385.     Icheck(r);
  1386.     rem = Itolong(r);
  1387.     delete r;
  1388.   }
  1389.   rem = abs(rem);
  1390.   if (xsgn == I_NEGATIVE) rem = -rem;
  1391.   q->sgn = samesign;
  1392.   Icheck(q);
  1393.   Iq.rep = q;
  1394. }
  1395.  
  1396.  
  1397. void divide(const Integer& Ix, const Integer& Iy, Integer& Iq, Integer& Ir)
  1398. {
  1399.   const IntRep* x = Ix.rep;
  1400.   nonnil(x);
  1401.   const IntRep* y = Iy.rep;
  1402.   nonnil(y);
  1403.   IntRep* q = Iq.rep;
  1404.   IntRep* r = Ir.rep;
  1405.  
  1406.   int xl = x->len;
  1407.   int yl = y->len;
  1408.   if (yl == 0)
  1409.     (*lib_error_handler)("Integer", "attempted division by zero");
  1410.  
  1411.   int comp = ucompare(x, y);
  1412.   int xsgn = x->sgn;
  1413.   int ysgn = y->sgn;
  1414.  
  1415.   int samesign = xsgn == ysgn;
  1416.  
  1417.   if (comp < 0)
  1418.   {
  1419.     q = Icopy_zero(q);
  1420.     r = Icopy(r, x);
  1421.   }
  1422.   else if (comp == 0)
  1423.   {
  1424.     q = Icopy_one(q, samesign);
  1425.     r = Icopy_zero(r);
  1426.   }
  1427.   else if (yl == 1)
  1428.   {
  1429.     q = Icopy(q, x);
  1430.     int rem =  unscale(q->s, q->len, y->s[0], q->s);
  1431.     r = Icopy_long(r, rem);
  1432.     if (rem != 0)
  1433.       r->sgn = xsgn;
  1434.   }
  1435.   else
  1436.   {
  1437.     IntRep* yy = 0;
  1438.     unsigned short prescale = (I_RADIX / (1 + y->s[yl - 1]));
  1439.     if (prescale != 1 || y == q || y == r)
  1440.     {
  1441.       yy = multiply(y, ((long)prescale & I_MAXNUM), yy);
  1442.       r = multiply(x, ((long)prescale & I_MAXNUM), r);
  1443.     }
  1444.     else
  1445.     {
  1446.       yy = (IntRep*)y;
  1447.       r = Icalloc(r, xl + 1);
  1448.       scpy(x->s, r->s, xl);
  1449.     }
  1450.  
  1451.     int ql = xl - yl + 1;
  1452.       
  1453.     q = Icalloc(q, ql);
  1454.     do_divide(r->s, yy->s, yl, q->s, ql);
  1455.  
  1456.     if (yy != y) delete yy;
  1457.     if (prescale != 1)
  1458.     {
  1459.       Icheck(r);
  1460.       unscale(r->s, r->len, prescale, r->s);
  1461.     }
  1462.   }
  1463.   q->sgn = samesign;
  1464.   Icheck(q);
  1465.   Iq.rep = q;
  1466.   Icheck(r);
  1467.   Ir.rep = r;
  1468. }
  1469.  
  1470. IntRep* mod(const IntRep* x, const IntRep* y, IntRep* r)
  1471. {
  1472.   nonnil(x);
  1473.   nonnil(y);
  1474.   int xl = x->len;
  1475.   int yl = y->len;
  1476.   if (yl == 0) (*lib_error_handler)("Integer", "attempted division by zero");
  1477.  
  1478.   int comp = ucompare(x, y);
  1479.   int xsgn = x->sgn;
  1480.  
  1481.   if (comp < 0)
  1482.     r = Icopy(r, x);
  1483.   else if (comp == 0)
  1484.     r = Icopy_zero(r);
  1485.   else if (yl == 1)
  1486.   {
  1487.     int rem = unscale(x->s, xl, y->s[0], 0);
  1488.     r = Icopy_long(r, rem);
  1489.     if (rem != 0)
  1490.       r->sgn = xsgn;
  1491.   }
  1492.   else
  1493.   {
  1494.     IntRep* yy = 0;
  1495.     unsigned short prescale = (I_RADIX / (1 + y->s[yl - 1]));
  1496.     if (prescale != 1 || y == r)
  1497.     {
  1498.       yy = multiply(y, ((long)prescale & I_MAXNUM), yy);
  1499.       r = multiply(x, ((long)prescale & I_MAXNUM), r);
  1500.     }
  1501.     else
  1502.     {
  1503.       yy = (IntRep*)y;
  1504.       r = Icalloc(r, xl + 1);
  1505.       scpy(x->s, r->s, xl);
  1506.     }
  1507.       
  1508.     do_divide(r->s, yy->s, yl, 0, xl - yl + 1);
  1509.  
  1510.     if (yy != y) delete yy;
  1511.  
  1512.     if (prescale != 1)
  1513.     {
  1514.       Icheck(r);
  1515.       unscale(r->s, r->len, prescale, r->s);
  1516.     }
  1517.   }
  1518.   Icheck(r);
  1519.   return r;
  1520. }
  1521.  
  1522. IntRep* mod(const IntRep* x, long y, IntRep* r)
  1523. {
  1524.   nonnil(x);
  1525.   int xl = x->len;
  1526.   if (y == 0) (*lib_error_handler)("Integer", "attempted division by zero");
  1527.  
  1528.   unsigned short ys[SHORT_PER_LONG];
  1529.   unsigned long u;
  1530.   int ysgn = y >= 0;
  1531.   if (ysgn)
  1532.     u = y;
  1533.   else
  1534.     u = -y;
  1535.   int yl = 0;
  1536.   while (u != 0)
  1537.   {
  1538.     ys[yl++] = extract(u);
  1539.     u = down(u);
  1540.   }
  1541.  
  1542.   int comp = xl - yl;
  1543.   if (comp == 0) comp = docmp(x->s, ys, xl);
  1544.  
  1545.   int xsgn = x->sgn;
  1546.  
  1547.   if (comp < 0)
  1548.     r = Icopy(r, x);
  1549.   else if (comp == 0)
  1550.     r = Icopy_zero(r);
  1551.   else if (yl == 1)
  1552.   {
  1553.     int rem = unscale(x->s, xl, ys[0], 0);
  1554.     r = Icopy_long(r, rem);
  1555.     if (rem != 0)
  1556.       r->sgn = xsgn;
  1557.   }
  1558.   else
  1559.   {
  1560.     unsigned short prescale = (I_RADIX / (1 + ys[yl - 1]));
  1561.     if (prescale != 1)
  1562.     {
  1563.       unsigned long prod = (unsigned long)prescale * (unsigned long)ys[0];
  1564.       ys[0] = extract(prod);
  1565.       prod = down(prod) + (unsigned long)prescale * (unsigned long)ys[1];
  1566.       ys[1] = extract(prod);
  1567.       r = multiply(x, ((long)prescale & I_MAXNUM), r);
  1568.     }
  1569.     else
  1570.     {
  1571.       r = Icalloc(r, xl + 1);
  1572.       scpy(x->s, r->s, xl);
  1573.     }
  1574.       
  1575.     do_divide(r->s, ys, yl, 0, xl - yl + 1);
  1576.  
  1577.     if (prescale != 1)
  1578.     {
  1579.       Icheck(r);
  1580.       unscale(r->s, r->len, prescale, r->s);
  1581.     }
  1582.   }
  1583.   Icheck(r);
  1584.   return r;
  1585. }
  1586.  
  1587. IntRep* lshift(const IntRep* x, long y, IntRep* r)
  1588. {
  1589.   nonnil(x);
  1590.   int xl = x->len;
  1591.   if (xl == 0 || y == 0)
  1592.   {
  1593.     r = Icopy(r, x);
  1594.     return r;
  1595.   }
  1596.  
  1597.   int xrsame = x == r;
  1598.   int rsgn = x->sgn;
  1599.  
  1600.   long ay = (y < 0)? -y : y;
  1601.   int bw = ay / I_SHIFT;
  1602.   int sw = ay % I_SHIFT;
  1603.  
  1604.   if (y > 0)
  1605.   {
  1606.     int rl = bw + xl + 1;
  1607.     if (xrsame)
  1608.       r = Iresize(r, rl);
  1609.     else
  1610.       r = Icalloc(r, rl);
  1611.  
  1612.     unsigned short* botr = r->s;
  1613.     unsigned short* rs = &(botr[rl - 1]);
  1614.     const unsigned short* botx = (xrsame)? botr : x->s;
  1615.     const unsigned short* xs = &(botx[xl - 1]);
  1616.     unsigned long a = 0;
  1617.     while (xs >= botx)
  1618.     {
  1619.       a = up(a) | ((unsigned long)(*xs--) << sw);
  1620.       *rs-- = extract(down(a));
  1621.     }
  1622.     *rs-- = extract(a);
  1623.     while (rs >= botr)
  1624.       *rs-- = 0;
  1625.   }
  1626.   else
  1627.   {
  1628.     int rl = xl - bw;
  1629.     if (rl < 0)
  1630.       r = Icopy_zero(r);
  1631.     else
  1632.     {
  1633.       if (xrsame)
  1634.         r = Iresize(r, rl);
  1635.       else
  1636.         r = Icalloc(r, rl);
  1637.       int rw = I_SHIFT - sw;
  1638.       unsigned short* rs = r->s;
  1639.       unsigned short* topr = &(rs[rl]);
  1640.       const unsigned short* botx = (xrsame)? rs : x->s;
  1641.       const unsigned short* xs =  &(botx[bw]);
  1642.       const unsigned short* topx = &(botx[xl]);
  1643.       unsigned long a = (unsigned long)(*xs++) >> sw;
  1644.       while (xs < topx)
  1645.       {
  1646.         a |= (unsigned long)(*xs++) << rw;
  1647.         *rs++ = extract(a);
  1648.         a = down(a);
  1649.       }
  1650.       *rs++ = extract(a);
  1651.       if (xrsame) topr = (unsigned short*)topx;
  1652.       while (rs < topr)
  1653.         *rs++ = 0;
  1654.     }
  1655.   }
  1656.   r->sgn = rsgn;
  1657.   Icheck(r);
  1658.   return r;
  1659. }
  1660.  
  1661. IntRep* lshift(const IntRep* x, const IntRep* yy, int negatey, IntRep* r)
  1662. {
  1663.   long y = Itolong(yy);
  1664.   if (negatey)
  1665.     y = -y;
  1666.  
  1667.   return lshift(x, y, r);
  1668. }
  1669.  
  1670. IntRep* bitop(const IntRep* x, const IntRep* y, IntRep* r, char op)
  1671. {
  1672.   nonnil(x);
  1673.   nonnil(y);
  1674.   int xl = x->len;
  1675.   int yl = y->len;
  1676.   int xsgn = x->sgn;
  1677.   int xrsame = x == r;
  1678.   int yrsame = y == r;
  1679.   if (xrsame || yrsame)
  1680.     r = Iresize(r, calc_len(xl, yl, 0));
  1681.   else
  1682.     r = Icalloc(r, calc_len(xl, yl, 0));
  1683.   r->sgn = xsgn;
  1684.   unsigned short* rs = r->s;
  1685.   unsigned short* topr = &(rs[r->len]);
  1686.   const unsigned short* as;
  1687.   const unsigned short* bs;
  1688.   const unsigned short* topb;
  1689.   if (xl >= yl)
  1690.   {
  1691.     as = (xrsame)? rs : x->s;
  1692.     bs = (yrsame)? rs : y->s;
  1693.     topb = &(bs[yl]);
  1694.   }
  1695.   else
  1696.   {
  1697.     bs = (xrsame)? rs : x->s;
  1698.     topb = &(bs[xl]);
  1699.     as = (yrsame)? rs : y->s;
  1700.   }
  1701.  
  1702.   switch (op)
  1703.   {
  1704.   case '&':
  1705.     while (bs < topb) *rs++ = *as++ & *bs++;
  1706.     while (rs < topr) *rs++ = 0;
  1707.     break;
  1708.   case '|':
  1709.     while (bs < topb) *rs++ = *as++ | *bs++;
  1710.     while (rs < topr) *rs++ = *as++;
  1711.     break;
  1712.   case '^':
  1713.     while (bs < topb) *rs++ = *as++ ^ *bs++;
  1714.     while (rs < topr) *rs++ = *as++;
  1715.     break;
  1716.   }
  1717.   Icheck(r);
  1718.   return r;
  1719. }
  1720.  
  1721. IntRep* bitop(const IntRep* x, long y, IntRep* r, char op)
  1722. {
  1723.   nonnil(x);
  1724.   unsigned short tmp[SHORT_PER_LONG];
  1725.   unsigned long u;
  1726.   int newsgn;
  1727.   if (newsgn = (y >= 0))
  1728.     u = y;
  1729.   else
  1730.     u = -y;
  1731.   
  1732.   int l = 0;
  1733.   while (u != 0)
  1734.   {
  1735.     tmp[l++] = extract(u);
  1736.     u = down(u);
  1737.   }
  1738.  
  1739.   int xl = x->len;
  1740.   int yl = l;
  1741.   int xsgn = x->sgn;
  1742.   int xrsame = x == r;
  1743.   if (xrsame)
  1744.     r = Iresize(r, calc_len(xl, yl, 0));
  1745.   else
  1746.     r = Icalloc(r, calc_len(xl, yl, 0));
  1747.   r->sgn = xsgn;
  1748.   unsigned short* rs = r->s;
  1749.   unsigned short* topr = &(rs[r->len]);
  1750.   const unsigned short* as;
  1751.   const unsigned short* bs;
  1752.   const unsigned short* topb;
  1753.   if (xl >= yl)
  1754.   {
  1755.     as = (xrsame)? rs : x->s;
  1756.     bs = tmp;
  1757.     topb = &(bs[yl]);
  1758.   }
  1759.   else
  1760.   {
  1761.     bs = (xrsame)? rs : x->s;
  1762.     topb = &(bs[xl]);
  1763.     as = tmp;
  1764.   }
  1765.  
  1766.   switch (op)
  1767.   {
  1768.   case '&':
  1769.     while (bs < topb) *rs++ = *as++ & *bs++;
  1770.     while (rs < topr) *rs++ = 0;
  1771.     break;
  1772.   case '|':
  1773.     while (bs < topb) *rs++ = *as++ | *bs++;
  1774.     while (rs < topr) *rs++ = *as++;
  1775.     break;
  1776.   case '^':
  1777.     while (bs < topb) *rs++ = *as++ ^ *bs++;
  1778.     while (rs < topr) *rs++ = *as++;
  1779.     break;
  1780.   }
  1781.   Icheck(r);
  1782.   return r;
  1783. }
  1784.  
  1785.  
  1786.  
  1787. IntRep*  compl(const IntRep* src, IntRep* r)
  1788. {
  1789.   nonnil(src);
  1790.   r = Icopy(r, src);
  1791.   unsigned short* s = r->s;
  1792.   unsigned short* top = &(s[r->len - 1]);
  1793.   while (s < top)
  1794.   {
  1795.     unsigned short cmp = ~(*s);
  1796.     *s++ = cmp;
  1797.   }
  1798.   unsigned short a = *s;
  1799.   unsigned short b = 0;
  1800.   while (a != 0)
  1801.   {
  1802.     b <<= 1;
  1803.     if (!(a & 1)) b |= 1;
  1804.     a >>= 1;
  1805.   }
  1806.   *s = b;
  1807.   Icheck(r);
  1808.   return r;
  1809. }
  1810.  
  1811. void setbit(Integer& x, long b)
  1812. {
  1813.   if (b >= 0)
  1814.   {
  1815.     int bw = (unsigned long)b / I_SHIFT;
  1816.     int sw = (unsigned long)b % I_SHIFT;
  1817.     if (x.rep == 0)
  1818.       x.rep = Icalloc(0, bw + 1);
  1819.     else if (x.rep->len < bw)
  1820.     {
  1821.       int xl = x.rep->len;
  1822.       x.rep = Iresize(x.rep, calc_len(xl, bw+1, 0));
  1823.     }
  1824.     x.rep->s[bw] |= (1 << sw);
  1825.     Icheck(x.rep);
  1826.   }
  1827. }
  1828.  
  1829. void clearbit(Integer& x, long b)
  1830. {
  1831.   if (b >= 0)
  1832.   {
  1833.     int bw = (unsigned long)b / I_SHIFT;
  1834.     int sw = (unsigned long)b % I_SHIFT;
  1835.     if (x.rep == 0)
  1836.       x.rep = Icalloc(0, bw + 1);
  1837.     else if (x.rep->len < bw)
  1838.     {
  1839.       int xl = x.rep->len;
  1840.       x.rep = Iresize(x.rep, calc_len(xl, bw+1, 0));
  1841.     }
  1842.     x.rep->s[bw] &= ~(1 << sw);
  1843.     Icheck(x.rep);
  1844.   }
  1845. }
  1846.  
  1847. int testbit(const Integer& x, long b)
  1848. {
  1849.   if (x.rep != 0 && b >= 0)
  1850.   {
  1851.     int bw = (unsigned long)b / I_SHIFT;
  1852.     int sw = (unsigned long)b % I_SHIFT;
  1853.     return (bw < x.rep->len && (x.rep->s[bw] & (1 << sw)) != 0);
  1854.   }
  1855.   else
  1856.     return 0;
  1857. }
  1858.  
  1859. // A  version of knuth's algorithm B / ex. 4.5.3.34
  1860. // A better version that doesn't bother shifting all of `t' forthcoming
  1861.  
  1862. IntRep* gcd(const IntRep* x, const IntRep* y)
  1863. {
  1864.   nonnil(x);
  1865.   nonnil(y);
  1866.   int ul = x->len;
  1867.   int vl = y->len;
  1868.   
  1869.   if (vl == 0)
  1870.     return Ialloc(0, x->s, ul, I_POSITIVE, ul);
  1871.   else if (ul == 0)
  1872.     return Ialloc(0, y->s, vl, I_POSITIVE, vl);
  1873.  
  1874.   IntRep* u = Ialloc(0, x->s, ul, I_POSITIVE, ul);
  1875.   IntRep* v = Ialloc(0, y->s, vl, I_POSITIVE, vl);
  1876.  
  1877. // find shift so that both not even
  1878.  
  1879.   long k = 0;
  1880.   int l = (ul <= vl)? ul : vl;
  1881.   int cont = 1;
  1882.   for (int i = 0;  i < l && cont; ++i)
  1883.   {
  1884.     unsigned long a =  (i < ul)? u->s[i] : 0;
  1885.     unsigned long b =  (i < vl)? v->s[i] : 0;
  1886.     for (int j = 0; j < I_SHIFT; ++j)
  1887.     {
  1888.       if ((a | b) & 1)
  1889.       {
  1890.         cont = 0;
  1891.         break;
  1892.       }
  1893.       else
  1894.       {
  1895.         ++k;
  1896.         a >>= 1;
  1897.         b >>= 1;
  1898.       }
  1899.     }
  1900.   }
  1901.  
  1902.   if (k != 0)
  1903.   {
  1904.     u = lshift(u, -k, u);
  1905.     v = lshift(v, -k, v);
  1906.   }
  1907.  
  1908.   IntRep* t;
  1909.   if (u->s[0] & 01)
  1910.     t = Ialloc(0, v->s, v->len, !v->sgn, v->len);
  1911.   else
  1912.     t = Ialloc(0, u->s, u->len, u->sgn, u->len);
  1913.  
  1914.   while (t->len != 0)
  1915.   {
  1916.     long s = 0;                 // shift t until odd
  1917.     cont = 1;
  1918.     int tl = t->len;
  1919.     for (i = 0; i < tl && cont; ++i)
  1920.     {
  1921.       unsigned long a = t->s[i];
  1922.       for (int j = 0; j < I_SHIFT; ++j)
  1923.       {
  1924.         if (a & 1)
  1925.         {
  1926.           cont = 0;
  1927.           break;
  1928.         }
  1929.         else
  1930.         {
  1931.           ++s;
  1932.           a >>= 1;
  1933.         }
  1934.       }
  1935.     }
  1936.  
  1937.     if (s != 0) t = lshift(t, -s, t);
  1938.  
  1939.     if (t->sgn == I_POSITIVE)
  1940.     {
  1941.       u = Icopy(u, t);
  1942.       t = add(t, 0, v, 1, t);
  1943.     }
  1944.     else
  1945.     {
  1946.       v = Ialloc(v, t->s, t->len, !t->sgn, t->len);
  1947.       t = add(t, 0, u, 0, t);
  1948.     }
  1949.   }
  1950.   delete t;
  1951.   delete v;
  1952.   if (k != 0) u = lshift(u, k, u);
  1953.   return u;
  1954. }
  1955.  
  1956.  
  1957.  
  1958. long lg(const IntRep* x)
  1959. {
  1960.   nonnil(x);
  1961.   int xl = x->len;
  1962.   if (xl == 0)
  1963.     return 0;
  1964.  
  1965.   long l = (xl - 1) * I_SHIFT - 1;
  1966.   unsigned short a = x->s[xl-1];
  1967.  
  1968.   while (a != 0)
  1969.   {
  1970.     a = a >> 1;
  1971.     ++l;
  1972.   }
  1973.   return l;
  1974. }
  1975.   
  1976. IntRep* power(const IntRep* x, long y, IntRep* r)
  1977. {
  1978.   nonnil(x);
  1979.   int sgn;
  1980.   if (x->sgn == I_POSITIVE || (!(y & 1)))
  1981.     sgn = I_POSITIVE;
  1982.   else
  1983.     sgn = I_NEGATIVE;
  1984.  
  1985.   int xl = x->len;
  1986.  
  1987.   if (y == 0 || (xl == 1 && x->s[0] == 1))
  1988.     r = Icopy_one(r, sgn);
  1989.   else if (xl == 0 || y < 0)
  1990.     r = Icopy_zero(r);
  1991.   else if (y == 1 || y == -1)
  1992.     r = Icopy(r, x);
  1993.   else
  1994.   {
  1995.     int maxsize = ((lg(x) + 1) * y) / I_SHIFT + 2;     // pre-allocate space
  1996.     IntRep* b = Ialloc(0, x->s, xl, I_POSITIVE, maxsize);
  1997.     b->len = xl;
  1998.     r = Icalloc(r, maxsize);
  1999.     r = Icopy_one(r, I_POSITIVE);
  2000.     for(;;)
  2001.     {
  2002.       if (y & 1)
  2003.         r = multiply(r, b, r);
  2004.       if ((y >>= 1) == 0)
  2005.         break;
  2006.       else
  2007.         b = multiply(b, b, b);
  2008.     }
  2009.     delete b;
  2010.   }
  2011.   r->sgn = sgn;
  2012.   Icheck(r);
  2013.   return r;
  2014. }
  2015.  
  2016. IntRep* abs(const IntRep* src, IntRep* dest)
  2017. {
  2018.   nonnil(src);
  2019.   if (src != dest)
  2020.     dest = Icopy(dest, src);
  2021.   dest->sgn = I_POSITIVE;
  2022.   return dest;
  2023. }
  2024.  
  2025. IntRep* negate(const IntRep* src, IntRep* dest)
  2026. {
  2027.   nonnil(src);
  2028.   if (src != dest)
  2029.     dest = Icopy(dest, src);
  2030.   if (dest->len != 0) 
  2031.     dest->sgn = !dest->sgn;
  2032.   return dest;
  2033. }
  2034.  
  2035. #if defined(__GNUG__) && !defined(NO_NRV)
  2036.  
  2037. Integer sqrt(const Integer& x) return r(x)
  2038. {
  2039.   int s = sign(x);
  2040.   if (s < 0) x.error("Attempted square root of negative Integer");
  2041.   if (s != 0)
  2042.   {
  2043.     r >>= (lg(x) / 2); // get close
  2044.     Integer q;
  2045.     div(x, r, q);
  2046.     while (q < r)
  2047.     {
  2048.       r += q;
  2049.       r >>= 1;
  2050.       div(x, r, q);
  2051.     }
  2052.   }
  2053.   return;
  2054. }
  2055.  
  2056. Integer lcm(const Integer& x, const Integer& y) return r
  2057. {
  2058.   if (!x.initialized() || !y.initialized())
  2059.     x.error("operation on uninitialized Integer");
  2060.   Integer g;
  2061.   if (sign(x) == 0 || sign(y) == 0)
  2062.     g = 1;
  2063.   else 
  2064.     g = gcd(x, y);
  2065.   div(x, g, r);
  2066.   mul(r, y, r);
  2067. }
  2068.  
  2069. #else 
  2070. Integer sqrt(const Integer& x) 
  2071. {
  2072.   Integer r(x);
  2073.   int s = sign(x);
  2074.   if (s < 0) x.error("Attempted square root of negative Integer");
  2075.   if (s != 0)
  2076.   {
  2077.     r >>= (lg(x) / 2); // get close
  2078.     Integer q;
  2079.     div(x, r, q);
  2080.     while (q < r)
  2081.     {
  2082.       r += q;
  2083.       r >>= 1;
  2084.       div(x, r, q);
  2085.     }
  2086.   }
  2087.   return r;
  2088. }
  2089.  
  2090. Integer lcm(const Integer& x, const Integer& y) 
  2091. {
  2092.   Integer r;
  2093.   if (!x.initialized() || !y.initialized())
  2094.     x.error("operation on uninitialized Integer");
  2095.   Integer g;
  2096.   if (sign(x) == 0 || sign(y) == 0)
  2097.     g = 1;
  2098.   else 
  2099.     g = gcd(x, y);
  2100.   div(x, g, r);
  2101.   mul(r, y, r);
  2102.   return r;
  2103. }
  2104.  
  2105. #endif
  2106.  
  2107.  
  2108.  
  2109. IntRep* atoIntRep(const char* s, int base)
  2110. {
  2111.   int sl = strlen(s);
  2112.   IntRep* r = Icalloc(0, sl * (lg(base) + 1) / I_SHIFT + 1);
  2113.   if (s != 0)
  2114.   {
  2115.     char sgn;
  2116.     while (isspace(*s)) ++s;
  2117.     if (*s == '-')
  2118.     {
  2119.       sgn = I_NEGATIVE;
  2120.       s++;
  2121.     }
  2122.     else if (*s == '+')
  2123.     {
  2124.       sgn = I_POSITIVE;
  2125.       s++;
  2126.     }
  2127.     else
  2128.       sgn = I_POSITIVE;
  2129.     for (;;)
  2130.     {
  2131.       long digit;
  2132.       if (*s >= '0' && *s <= '9') digit = *s - '0';
  2133.       else if (*s >= 'a' && *s <= 'z') digit = *s - 'a' + 10;
  2134.       else if (*s >= 'A' && *s <= 'Z') digit = *s - 'A' + 10;
  2135.       else break;
  2136.       if (digit >= base) break;
  2137.       r = multiply(r, base, r);
  2138.       r = add(r, 0, digit, r);
  2139.       ++s;
  2140.     }
  2141.     r->sgn = sgn;
  2142.   }
  2143.   return r;
  2144. }
  2145.  
  2146.  
  2147.  
  2148. extern AllocRing _libgxx_fmtq;
  2149.  
  2150. char* Itoa(const IntRep* x, int base, int width)
  2151. {
  2152.   int fmtlen = (x->len + 1) * I_SHIFT / lg(base) + 4 + width;
  2153.   char* fmtbase = (char *) _libgxx_fmtq.alloc(fmtlen);
  2154.   char* f = cvtItoa(x, fmtbase, fmtlen, base, 0, width, 0, ' ', 'X', 0);
  2155.   return f;
  2156. }
  2157.  
  2158. ostream& operator << (ostream& s, const Integer& y)
  2159. {
  2160. #ifdef _OLD_STREAMS
  2161.   return s << Itoa(y.rep);
  2162. #else
  2163.   int base = (s.flags() & ios::oct) ? 8 : (s.flags() & ios::hex) ? 16 : 10;
  2164.   int align_right = !(s.flags() & ios::left);
  2165.   int showpos = s.flags() & ios::showpos;
  2166.   int showbase = s.flags() & ios::showbase;
  2167.   char fillchar = s.fill();
  2168.   char Xcase = (s.flags() & ios::uppercase)? 'X' : 'x';
  2169.   int width = s.width();
  2170.   const IntRep* x = y.rep;
  2171.   int fmtlen = (x->len + 1) * I_SHIFT / lg(base) + 4 + width;
  2172.   char* fmtbase = new char[fmtlen];
  2173.   char* f = cvtItoa(x, fmtbase, fmtlen, base, showbase, width, align_right,
  2174.                     fillchar, Xcase, showpos);
  2175.   s.write(f, fmtlen);
  2176.   delete fmtbase;
  2177.   return s;
  2178. #endif
  2179. }
  2180.   
  2181.  
  2182. char*  cvtItoa(const IntRep* x, char* fmt, int& fmtlen, int base, int showbase,
  2183.               int width, int align_right, char fillchar, char Xcase, 
  2184.               int showpos)
  2185. {
  2186.   char* e = fmt + fmtlen - 1;
  2187.   char* s = e;
  2188.   *--s = 0;
  2189.  
  2190.   if (x->len == 0)
  2191.     *--s = '0';
  2192.   else
  2193.   {
  2194.     IntRep* z = Icopy(0, x);
  2195.  
  2196.     // split division by base into two parts: 
  2197.     // first divide by biggest power of base that fits in an unsigned short,
  2198.     // then use straight signed div/mods from there. 
  2199.  
  2200.     // find power
  2201.     int bpower = 1;
  2202.     unsigned short b = base;
  2203.     unsigned short maxb = I_MAXNUM / base;
  2204.     while (b < maxb)
  2205.     {
  2206.       b *= base;
  2207.       ++bpower;
  2208.     }
  2209.     for(;;)
  2210.     {
  2211.       int rem = unscale(z->s, z->len, b, z->s);
  2212.       Icheck(z);
  2213.       if (z->len == 0)
  2214.       {
  2215.         while (rem != 0)
  2216.         {
  2217.           char ch = rem % base;
  2218.           rem /= base;
  2219.           if (ch >= 10)
  2220.             ch += 'a' - 10;
  2221.           else
  2222.             ch += '0';
  2223.           *--s = ch;
  2224.         }
  2225.         delete z;
  2226.         break;
  2227.       }
  2228.       else
  2229.       {
  2230.         for (int i = 0; i < bpower; ++i)
  2231.         {
  2232.           char ch = rem % base;
  2233.           rem /= base;
  2234.           if (ch >= 10)
  2235.             ch += 'a' - 10;
  2236.           else
  2237.             ch += '0';
  2238.           *--s = ch;
  2239.         }
  2240.       }
  2241.     }
  2242.   }
  2243.  
  2244.   if (base == 8 && showbase) 
  2245.     *--s = '0';
  2246.   else if (base == 16 && showbase) 
  2247.   { 
  2248.     *--s = Xcase; 
  2249.     *--s = '0'; 
  2250.   }
  2251.   if (x->sgn == I_NEGATIVE) *--s = '-';
  2252.   else if (showpos) *--s = '+';
  2253.   int w = e - s - 1;
  2254.   if (!align_right || w >= width)
  2255.   {
  2256.     while (w++ < width)
  2257.       *--s = fillchar;
  2258.     fmtlen = e - s - 1;
  2259.     return s;
  2260.   }
  2261.   else
  2262.   {
  2263.     char* p = fmt;
  2264.     int gap = s - p;
  2265.     for (char* t = s; *t != 0; ++t, ++p) *p = *t;
  2266.     while (w++ < width) *p++ = fillchar;
  2267.     *p = 0;
  2268.     fmtlen = p - fmt;
  2269.     return fmt;
  2270.   }
  2271. }
  2272.  
  2273. char* dec(const Integer& x, int width)
  2274. {
  2275.   return Itoa(x, 10, width);
  2276. }
  2277.  
  2278. char* oct(const Integer& x, int width)
  2279. {
  2280.   return Itoa(x, 8, width);
  2281. }
  2282.  
  2283. char* hex(const Integer& x, int width)
  2284. {
  2285.   return Itoa(x, 16, width);
  2286. }
  2287.  
  2288. istream& operator >> (istream& s, Integer& y)
  2289. {
  2290. #ifdef _OLD_STREAMS
  2291.   if (!s.good())
  2292.     return s;
  2293. #else
  2294.   if (!s.ipfx(0))
  2295.   {
  2296.     s.set(ios::failbit);
  2297.     return s;
  2298.   }
  2299. #endif
  2300.   s >> ws;
  2301.   if (!s.good())
  2302.   {
  2303.     s.set(_fail);
  2304.     return s;
  2305.   }
  2306.   
  2307. #ifdef _OLD_STREAMS
  2308.   int know_base = 0;
  2309.   int base = 10;
  2310. #else
  2311.   int know_base = s.flags() & (ios::oct | ios::hex | ios::dec);
  2312.   int base = (s.flags() & ios::oct) ? 8 : (s.flags() & ios::hex) ? 16 : 10;
  2313. #endif
  2314.  
  2315.   int got_one = 0;
  2316.   char sgn = 0;
  2317.   char ch;
  2318.   y.rep = Icopy_zero(y.rep);
  2319.  
  2320.   while (s.get(ch))
  2321.   {
  2322.     if (ch == '-')
  2323.     {
  2324.       if (sgn == 0)
  2325.         sgn = '-';
  2326.       else
  2327.         break;
  2328.     }
  2329.     else if (!know_base & !got_one && ch == '0')
  2330.       base = 8;
  2331.     else if (!know_base & !got_one && base == 8 && (ch == 'X' || ch == 'x'))
  2332.       base = 16;
  2333.     else if (base == 8)
  2334.     {
  2335.       if (ch >= '0' && ch <= '7')
  2336.       {
  2337.         long digit = ch - '0';
  2338.         y <<= 3;
  2339.         y += digit;
  2340.         got_one = 1;
  2341.       }
  2342.       else
  2343.         break;
  2344.     }
  2345.     else if (base == 16)
  2346.     {
  2347.       long digit;
  2348.       if (ch >= '0' && ch <= '9')
  2349.         digit = ch - '0';
  2350.       else if (ch >= 'A' && ch <= 'F')
  2351.         digit = ch - 'A' + 10;
  2352.       else if (ch >= 'a' && ch <= 'f')
  2353.         digit = ch - 'a' + 10;
  2354.       else
  2355.         digit = base;
  2356.       if (digit < base)
  2357.       {
  2358.         y <<= 4;
  2359.         y += digit;
  2360.         got_one = 1;
  2361.       }
  2362.       else
  2363.         break;
  2364.     }
  2365.     else if (base == 10)
  2366.     {
  2367.       if (ch >= '0' && ch <= '9')
  2368.       {
  2369.         long digit = ch - '0';
  2370.         y *= 10;
  2371.         y += digit;
  2372.         got_one = 1;
  2373.       }
  2374.       else
  2375.         break;
  2376.     }
  2377.     else
  2378.       abort(); // can't happen for now
  2379.   }
  2380.   if (s.good())
  2381.     s.putback(ch);
  2382.   if (!got_one)
  2383.     s.set(_fail);
  2384.  
  2385.   if (sgn == '-')
  2386.     y.negate();
  2387.  
  2388.   return s;
  2389. }
  2390.  
  2391. int Integer::OK() const
  2392. {
  2393.   int v = rep != 0;             // have a rep
  2394.   int l = rep->len;
  2395.   int s = rep->sgn;
  2396.   v &= l <= rep->sz;            // length with bounds
  2397.   v &= s == 0 || s == 1;        // legal sign
  2398.   Icheck(rep);                  // and correctly adjusted
  2399.   v &= rep->len == l;
  2400.   v &= rep->sgn == s;
  2401.   if (!v) error("invariant failure");
  2402.   return v;
  2403. }
  2404.  
  2405. void Integer::error(const char* msg) const
  2406. {
  2407.   (*lib_error_handler)("Integer", msg);
  2408. }
  2409.  
  2410.